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Abstract 


New technological needs are leading to the development of new materials with enhanced 
mechanical properties for specific applications. The knowledge of the elastic properties 
of these new materials, i.e. the elastic constants, is in practice of great importance. Many 
of these new materials are ceramics and composites. Conventional methods are not very 
much successful for determining these material properties. 

The laser ultrasonic technique, which uses lasers to generate and detect ultrasound at a 
distance, without any contact to the specimen, are more popular and successful for such 
materials. If the samples can be cut in predetermined directions, elastic constants can be 
determined from velocity data in these directions with the number of measurements equal 
to number of unknown independent elastic constants. Choosing an appropriate 
propagation direction can maximize the accuracy of the elastic constant measurement. 
The process of image reconstruction using velocity data is termed as tomography. 

In the present work, a tomographic image reconstruction model is developed using ray 
tracing. The method reconstructs a two dimensional object, with internal defects, in terms 
of independent parameters, like elastic constants of the material. Single mode 
measurement time of flight data is used for the reconstruction. The proposed model is a 
two stage process, which involves conventional ART algorithm corrections in the first 
stage and solution of a set of nonlinear equations in the second stage. For defects like 
circle and ellipse, edge detection is taken into consideration, during ray tracing. The 
method is found to perform well for rectangular, as well as curved defects. ' 
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Chapter 1 

Introduction 


1.1 INTRODUCTION 

The search for ever lighter, stiffer, and tougher materials from which to construct aerospace 
vehicles has led inexorably towards more highly engineered, lower density materials like 
cornposite materials. The novelty and expense of composite materials has meant, historically, 
that they are thoroughly inspected prior to use, typically using ultrasonic non-destructive 
testing. 

Ultrasonic materials characterization is the most important application of ultrasonics in 
aerospace engineering and engineering mechanics. Historically, ultrasonic nondestructive 
testing (NDT) has been used almost exclusively for detecting macroscopic discontinuities in 
structures after they have been in service for some time. It has become increasingly evident 
that it is practical and cost effective to expand the role of ultrasonic NDT testing to include 
all aspects of materials production and application. Research efforts are being directed at 
developing and perfecting NDT capable of monitoring (i) material production processes, (ii) 
material integrity following transport, storage and fabrication, and (iii) the amount and rate of 
degradation during service. In addition, efforts are underway to develop techniques capable 
of quantitative discontinuity sizing, permitting determination of material response using 
fracture mechanics analysis, as well as techniques for quantitative materials characterization 
to replace the qualitative techniques used in the past. Ultrasonic techniques play a prominent 
role in these developments because they afford useful and versatile methods for evaluating 
microstructures, associated mechanical properties, as well as detecting microscopic and 
macroscopic discontinuities in solid materials. 

Ultrasonic measurements are widely used for elastic constant determination of anisotropic 
materials. Tomography is a powerful NDE tool that has seen a rapid development in last two 
decades. Tomography is image reconstruction of functions from multiple projections. In 
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tomography, two basic approaches are used to obtain the information required to reconstruct 
the generally complex propagation velocity distribution. These approaches can be classified 
as ray-tracing methods and wave-propagation methods. Both methods are based on 
approximate solutions of the wave equation, which governs the propagation of ultrasonic 
wave through media with inhomogeneous velocity distribution. 

In Tomographic reconstruction, an object having some property distribution is inspected 
from a number of angles. These projections are then used to generate a distribution over a 
grid of cells, referred to as pixels, by superposition of all the projection at each pixel. 
Basically there are two types of reconstruction methods, viz., direct and iterative. The direct 
method can be further classified into direct algebraic methods and Fourier transform based 
methods. The iterative methods include various reconstruction techniques such as Algebraic 
Reconstruction Technique (ART), Simultaneous Iterative Reconstruction Technique (SIRT), 
Multiplicative Algebraic Reconstruction Technique (MART), and Entropy Optimization etc. 
These reconstruction algorithms are frequently applied to determine elastic constants of 
composite materials. 

1.2 LITERATURE SURVEY 

The conventional tomographic problem may be described as the task of reconstructing an 
unknown scalar function from its projections. In a variety of applications, it has been 
successfully applied to get a scalar property distribution in a cross section, such as slowness 
distribution etc. 

Synge [1] developed the theory of slowness surfaces after that Musgrave [2] first introduced 
the idea of superimposing the slowness surface on the geometric space, which contains the 
incident and emerging waves. 

In the literature, only a few ray-tracing algorithm for wave propagation in anisotropic 
inhomogeneous materials are discussed. Silk [3] divides the austenitic weld region into a 
number of quadrilateral regions, each with a fixed direction of grain crystals. Ogilvy [4] 
defined the inhomogeneity of the weld by empirical analytic function that describes the 
orientation of grain crystals continuously. 
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More recently, Rokhlin et al. [5] have described a unified approach to the numerical solution 
of the reflection-refraction problem for generally anisotropic media. Ogilvy [6] also studied 
the influence of austenitic weld geometry and manufacture on ultrasonic inspection of weld 
joints. Ogilvy [7,8] was the first to purpose a ray-tracing model for ultrasonic NDT in 
austenite welds. Schmitz et al. [9] modified the Ogilvy’s model and proposed a new model 
for austenite materials to follow the longitudinal, horizontal and vertical polarized shear 
wave propagation from the base material through the cladding in three dimensions. 

The treatment of anisotropy has been dealt with, by some researchers, for various 
applications. Thomson [10] presented simplified expressions for wave velocities in terms of 
anisotropic parameters and trigonometric functions. These expressions are appropriate for 
weakly anisotropic media (10-20 percent), occurring in geophysics applications. This weak 
anisotropy formulation was used by Stewart [11] in deriving Algebraic Reconstruction 
Technique (ART) correction for base velocity and anisotropic parameters. These ART 
corrections are applicable in case of cross hole imaging. This approach was later applied for 
tomographic reconstruction in rocks and polymer composites [12,13]. Kline and Wang [14] 
developed a modified versions of ART for composite materials. They used the governing 
equations for wave propagation in anisotropic material and composite micro mechanism to 
account for variation of acoustic velocity with propagation direction. This method requires a 
prior knowledge of properties of composite material constituents, that is, fibers and matrix. 
The time delay measurement of single wave mode was used for tomographic reconstruction. 
This technique was later extended to construct multiple material parameters using the time 
delay information of multiple independent modes of acoustic signal [15]. Generally when the 
time delay measurement multiple modes are considered, it is required to perform ray tracing 
for each mode independently. This makes it computationally intentional, as the method is 
iterative. 

1.3 PRESENT WORK 

In this work, ray tracing technique is applied for tomographic image reconstruction to 
determine principal elastic constants and density of given specimen. Multi-parameter image 
reconstruction using ART algorithm is presented. 
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To ensure the reliability of ultrasonic inspection techniques, these complex material 
properties as well as the influence of micro structural inhomogeneities and the effects of non- 
planar surfaces and interfaces on ultrasonic wave propagation have to be taken in to account. 
In this respect, simulation and optimization in ultrasonic testing have gained a considerable 
importance, where mathematical modeling provides a simple method of assisting analysis. In 
this presentation a ray-tracing model is developed to follow quasi-longitudinal and quasi- 
shear wave propagation in composite materials for various phase angle, which is further 
applied to determine elastic constants and density of composite materials. 

The present model is also used to reconstruct images using data obtained from Laser Based 
Ultrasonic (LBU) setup. Effect of various parameters on the performance of reconstruction 
model is evaluated. 

1.4 THESIS ORGANISATION 

Chapter II gives the fundamentals of elastic wave propagation in anisotropic media and 
refraction principles based on Snell’s law. It describes the equations of motion and plane 
wave solution of this equation. It discusses the fundamental equation of wave propagation in 
inhomogeneous and anisotropic media. Various terms associated with the elastic wave 
propagation in anisotropic media such as phase and group velocity, slowness, energy flux 
and critical angle are discussed in this chapter. Snell’s law applied to anisotropic materials is 
also discussed. 

Chapter III presents mathematical formulation of the ray tracing and then algorithm of two 
point ray tracing is discussed in detail. Then application of Edge detection on ray tracing 
algorithm is discussed. 

Chapter IV presents the basic tomographic algorithms. Then application of tomography on 
determination of elastic constants and density is discussed. 

Chapter V includes the results of application of tomographic image reconstruction model in 
determining elastic constants. Then experimental results based on experimental data obtained 
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from Laser ultrasonic set up are discussed. Then effect of various parameters on the 
performance of model is discussed. 

Chapter VI presents the conclusion and scope of future work. 
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Chapter 2 

Basics of Ultrasonic wave propagation in anisotropic 

media 


2.1 INTRODUCTION 

Ultrasonics involves the propagation of acoustic waves. Therefore, it is necessary to 
understand the basic features of propagating waves and some of the mathematical equations 
governing wave propagation. Acoustics is the study of time-varying deformations, or 
vibrations in elastic media. It is concerned with material particles that are small but yet 
contain many atoms. Within each particle the atoms move in unison. Therefore, acoustics 
deals with macroscopic phenomena and is formulated as if matter were a continuum. 

Propagation of ultrasonic waves in isotropic media has been well studied and can be 
predicted theoretically with good results. The problem is defined by two constants for the 
medium, such as X and p, the Lame constants and, the incident wave speed and direction. 
For anisotropic media the problem becomes considerably more complicated, as the elastic 
constants determining the wave propagation are then directionally dependent. 

In general, for a given propagation direction, the wave equation has three solutions with 
dilTcring speeds corresponding to three different polarization. A plane wave solution to the 
wave equation describes the possibility of three different wave modes with different phase 
velocities, 't hese wave modes are not in general, purely longitudinal or transverse, but quasi 
wave modes. 

The effect of anisotropy on ultrasotiic wave propagation can be very marked as wave 
velocities become directionally dependent; group and phase velocities are no longer 
necessarily parallel or equal in magnitude, 'fhese effects lead to beam distortion with in the 
material, unexpected beam propagation directions, ‘favoured’ propagation directions, 
alteration in spatial profiles and, for inhomogeneous media, beam paths which become 
curved as the beam propagates through the material. It is therefore important to be able to 
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predict such effects to ensure that beam energy (wliich travels in the group velocity 
direction) is directed to tire required region for inspection purpose. A brief review of wave 
propagation is presented. The basic concepts and formulations, included in this chapter are 
taken from[17J. 

2.2 KQDA I ION OF MOTION 

The static equilibrium of an arbitrary, three-dimensional solid body occupying volume V 
that has both surface and body forces requires that the resultant force vector vanish, this 
gives 


‘^ + f?b,=Q /=l,2,3 (2.1) 

dXj 

where Oij represents components of stress and p is density and bi is the body force density 
vector.'! his is a static result, however, and to describe the dynamical motion of a continuum, 
an inertial term accounting for the particle acceleration must be included. 


da 3 = 1 / 

‘+pb, = p- 


dx, 


dV 


( 2 . 2 ) 


It is also called Cauchy's equation. If body force hi vanishes for free, unforced motion, it 
gives 


Od,. d^u. 

3.; 


(2,3) 


and this gives dynamical equation of equilibrium for the elastic continuum. Using 
constitutive relations Hq. 2.3 can be written as 


r ^ 
dx, 


d\ 

de 


(2.4) 


where, Cijki is the elastic coefficient tensor. It is symmetric and positive definite in the sense 
that 
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for all symmetric tiij where equality is satisfied only when nij= 0 

By inserting for the strain Cki its linearized equivalent in terms of the displacement gradient, 


Eq. (2.4) can be written as 


1 his can be written as 


dx, 


■ + 


du, 

dx. 


1 ^ d du, ^ du, 




dx,dx^ 


= P 


d\ 

dP 


(2.5) 


( 2 . 6 ) 


(2.7) 


Equation (2.7) is the equation of motion for an infinite linear elastic medium in three 
dimensions. 


2.3 PLANE WAVE PROPAGATION IN BULK MATERIALS 


The plane harmonic wave solution of Eq. (2.7) can be expressed either in index notation or 
in vector notation as 


u(r,/) = /l„pexp/(k.r- tyf) 

(2.8) 

- ax) 

(2.9) 


where Ao is the amplitude of the wave, p (pj) is the polarization unit vector and co is the 
angular frequency. The wavevector k or (km) is given by 


k = A:ii, 



( 2 . 10 ) 


where n is a unit vector in the direction of wave propagation, and X is the wavelength. The 
phase of the wave (j) is given by 
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= k ■ r - liV 


( 2 . 11 ) 


For points of constant phase on the wave is 

\i.-r - ojt = comt. (2.12) 

Taking differential and rearranging gives 


k- — = (2.13) 

di 

Substituting k from Eq. (2. 10) and let ^ be the projection of the wavefront position vector r 
on the direction of wave propagation n . r), as shown in Fig. 2.1, 


All • 


di 


= 0 ) 


(2.14) 





Figure 2. 1 Wavefiont moving in a rectangular coordinate system 


and the propagation of a point of constant phase proceeds with pliase velocity V given by 


di k 


(2.15) 


in vectorial form 
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(2.16) 


V = rii = -ii. 
k 

Subslituling the plane harmonic wave solution (2.9) into the equation of motion (2.7) gives ^ 

(2.17) 


Or 




(2.18) 


By using above equation can be rewritten as 




(2.19) 


where 6ini is the Kroneker delta. In order to have a nontrivial solution for pm, the detenninant 
of the coefficients for the linear equation Eq. (2.18) must vanish, leading to the following 
well known eigenvalue equation, also known as Christoirel's equation. 




( 2 . 20 ) 


if A,ii = Ciji,„njni and ti = pV^ 

Hq. (2. 1 9) can also be written as 


(>^ii - B^im) Pm=0 (2.21) 

From the properties of Cijim, it follows that A-n is also symmetric and positive definite. 
Therefore, all the Eigen values of A,ii are real and positive and their corresponding 
eigenvectors are orthogonal. 
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This cat! he in(ci[)ie(c(J as that a Cor a given diicclioii of wave piopugatioii n there will be 
three phase velocities, and three corresponding displacement vectors will be orthogonal. 
These three waves are pure shear, quasi-longitudinal and quasi shear wave. 


Pure longitudinal mode (the polarization vector p is parallel to the propagation direction m, 
1, n are the direction cosines of propagation direction) is obtained only when p,=n, (i = 
1,2,3). Pure shear mode (the polarization vector p is normal to the propagation direction n) 
are obtained only when |> n = 0) In all other cases, either quasi- 

longitudinal or quasi-shear modes can only propagate. It can also be noted that any general 
wave propagation can be split in three modes: pure shear (S), quasi-longitudinal (QL) and 
quasi-shcar (QS). In a elaborate form Eq.( 2.21) can be written as 



^.2 

^.3 

'P^' 


0 

^2 

A 22 - pV^ 

^^2.3 

' P 2 

^ 

0 

^.3 

^22 

A,,-pV\ 

.AJ 


0 


( 2 . 22 ) 


where X in terms of elastic constants can be given as 
;ii,=^^C,,,,+w^C,2,2+^»^C',„2+2/»/2C,3,2+2«fC,,„+2mfC,,,2 

^22=^' +2mnC,,,, -f 2/2t'C2,,2 +2wfC22,2 

f ^ C,„ , + 2/»22C„2., +2n(:C„„ +2/;;fC2„ , (2.23) 

— ^ C,,, 2+222 ^2212 +22 C'23|3 +22222(C23|2+C22i3) + 22f(C,|23+C,j|2)+222£(C|,22+C'j2|2) 

= ( C, , , 3 + 222 C 23 1 2 + 2 2 C 33 1 3 + 22 222(^23 1 3 + 1 2 )+ 22 f ((7, , 33 + C", 3 , 3 ) + 222 ^((7, , 23 + C| 3 1 2 ) 

^23 = £ C,3,2 +222 C2223 +22 ^^3323 2.32.3 '*'^22.3.3 .3.312 '^'^'2.3 1.3)’^ ^^^^(^22 1 3 '^^2312) 

m, I, n are the direction cosines of propagation direction 


2.4 PHASE AND GROUP VELOCITY 


Pha.se velocity is defined as the velocity with which plane wave crests and troughs travel 
through a medium and is expressed as the ratio of the frequency of vibration and the wave 
number (i.e. the number of wavelengths per unit distance normal to the wavefronts). Group 


II 



velocity, also known as ray or energy velocity, is delincd as the velocity with which the 
energy of the wave propagates. 

An additional complication in anisotropic media is that the wave normal does not 
necessarily coincide with the propagation direction of the energy (ray direction)^ An 
arbitrary wave is dependent on time and spatial coordinates as / = f{G3(-kx}. In 
dispersive media, where the velocity is dependent on frequency, the acoustic energy of a 
narrow-band tone-burst propagates at the group velocity Vg. The phase velocity is simply 


and the group velocity is given as 



V = — 
* dk 


(2.24) 


(2.25) 


Transmitter 


Beam Contour 


Phase 



Figure 2.2 Deviaiton between phase direction and ray direction 


In the symmetry directions only pure modes propagate and the phase and ray directions are 
necessarily parallel. I his is because deviation from a symmetry direction in any direction 
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itiusl Imve the same cllcct as an equal deviation in llic opposite direction. This also means 
that the velocities of pure waves are either maxima or minima for that particular mode. In 
any other directions the beam becomes skewed as illustrated in Figure 2.2. Needless to say 
that this skewing effect makes it more difficult to align the ultrasonic transducers for optimal 
inspection in highly anisotropic materials such as unidirectional composite laminates. The 
same effect is also responsible for the lateral shift of the transmitted ultrasonic beam through 
anisotropic plates that occurs even at normal incidence. 

2.5 RELATION BETWEEN THE PHASE AND GROUP VELOCITIES 

The group velocity can be written in terms of phase velocity as 

" pV 

In order to find the relation between Vg and V a vector Q is defined here as follows 

Q,=C,„„n,p„p, (2.27) 


(2.26) 


group velocity can be written as 


(2.28) 

pV 


From the equation of motion 


Cy,„njn,p„p,-pV^=^^. (2.29) 

This can be written as 

Q,n, = pV^ (2.30) 

Finally, from Eqs. (2.46) and (2.48), 


/ 
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Vgcosy^y, 


(2.32) 


where y is the angle between and V. I his relation indicates that V is the projection of Vg 
in the direction of wave propagation 



Figure 2.3 Relation beween Phase velocity vector (jroi^ velocity vector 


2.6 SLOWNESS SURFACE 


An extremely useful insight of the effects of anisotropy can be obtained by examination of 
the slowness surfaces of a material. The correct interpretation of these surfaces can aid in the 
prediction of beam profiles in anisotropic materials. 


The slowness vector of a wave, ni, is defined as 


Ilciicc 111 hns mngniludc given by Ihc reciprocal of the phase velocity. Since, for an 
anisotropic material, a wave’s phase velocity magnitude is dependent on the phase velocity 
direction, it can be said that |m| is a function of n, i.e. 

Ill = m(n) (2.34) 

If the value of ni is plotted for all possible values of n, i.e. all phase velocity direction, then 
a three-dimensional surface mapped out by the end points of the slowness vectors, is called 
slowness surface. Since there are three modes of wave propagation in anisotropic materials 
there exists three slowness surfaces corresponding to the three wave modes. 

For an isotropic material each surface, corresponding to a wave polarization, will be a 
sphere since m is independent of n (i.e. phase velocity is independent of direction in 
isotropic materials) 

Figure 2.6 shows a slice through the slowness surface for transversely isotropic. There are 
several points to notice from the figure 

1. Each slowness surface is not a circle, indicating departure from isotropy. 

2. The three slowness surfaces, one for each mode, are of difierent shapes. This immediately 
indicates that no two wave modes will behave in the same manner with in the composite. 

A further use of slowness surface arises in depicting group velocity i.e. energy direction. 
Theory shows that for any phase velocity direction n, the normal to slowness surface for that 
particular value of ii is in the direction of the group velocity. 

The usefulness of this unusual parameter than conventional velocity, lies in the nature of 
anisotropy itself. In many cases (e. g., refraction at a plane interface or diffraction from a 
periodic structure) the wave direction is determined by the wave speed. For example, 
according to the well-known Snell's law, the refraction angles of all reflected and 
transmitted waves are determined by the incident angle 0/ and incident velocity as 
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follows 


Figure 2.4 A general slowness diagram for an anisotropic material 



' (2.35) 

n K 

where the subscript r represents any of the refracted waves. This simple explicit relation for 
0^ becomes an implicit one for anisotropic cases 

sin^^ _sin^, 

With slowness instead of velocity, the above condition becomes 

{0^ )sin 0^ = m, {0, )sin 0, , 

i.e., the projections of the slownesses of all waves are the same on the interface. This feature 
can be readily exploited to determine the refraction angles in a graphical way. At first the 
horizontal slowness projection nt = m, {0, )sin 0, is determined from the angle of incidence. 


(2.36) 


(2.37) 

i 
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Then, for ench particular rclYaclcd wave the refraction angle is determined from the 
intersection point of the slowness diagram with the same vertical line representing m. A 
special case in which the vertical line does not intersect a particular mode’s slowness 
diagram then that mode will no longer exist in propagation. 



Figure 2.5 representation ofreflected/transnutted wave modes using slowness diagram 

2.7 FUNDAMEN l ALS OF REFLECTION AND REFRACTION 

The simplest situation is depicted in Figure 2.6a, where a wave encounters a boundary at 
right angle or normal incidence. The interaction only involves reflection of some of the 
wave and transmission of a portion, with the amount of energy in each part depending on the 
material characteristics. A more complicated situation may arise, particularly in solids, when 
the wave strikes at an angle, or at oblique incidence. What may occur, as shown in Fig. 2.6b, 
is that two types of waves are reflected for a single incident wave. This phenomenon is 
known as mode conversion, and is illustrated for the case of a longitudinal wave generkting 
both longitudinal and shear waves. Yet another aspect is involved when waves encounter 
edges. Complex scattering and diffraction of the waves may occur, similar to optics (figure 
2.6c). 
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Figure 2.6 DiflEerent brpe of acoustic wave interaction with material 
discontinuties; (a) Normal Incidence (b) Oblique Incidence and Mode 
conversion (c) Diffiaction and Scattermg 
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2.8 SNICLL’S LAW 


Application of Snell’s law for anisotropic materials can be understand by considering a 
monochromatic plane wave incident on the interface between two anisotropic solids as 
shown in Fig. 2.7 the interface coincide with X 1 X 2 plane and b is a unit vector parallel to the 
X3. 

Tlie displacement vector uj of an elastic wave can be given as 1 

Uj = A^Pj exp[/(A:,x, - ox)] (2.38) 

In tenns of the previously defined slowness vector m, the displacement vector can be also 
written in the following alternative form 

Uj = APj exp[/ry(m,JC, - 1)] (2.39) 



Figure 2 . 1 Acoustic wave incident on aji interface betweeen two 
anisotropic media 


Because of assumed rigid boundary conditions at the interface, it requires continuity of both 
displacements and tractiotis across the boundary. Then, continuity of the displacements and 
tractions at the boundary (X 3 “ 0 ) is given by 


19 



and 


tr-l n-A 

(2.40) 

«-! a^A 

(2.41) 


1 he incident wave is marked with superscript 0 , and the three reflected and three transmitted 
waves are marked with superscripts 1 ; 2 ; 3 and 4 ; 5 ; 6, respectively. 


Using the relevant displacement-strain and constitutive relationships, the stress tensor can be 
expressed as 


du 


Qx\{io){m,x, - /)J (2.42) 


dx 


Substituting [!q. (2.39) into (2.40) the displacement continuity (2.40) can be written as ) 
Kp] exp(/V^;?X)= exp{ia)itj;x;)+ ZA"pJ exp(/aMn;‘x"), (2.43) 


Similarly Substituting Eq. (2.42) into (2.41) the stress continuity (2.41) can be written as 
exp(/rMr??;jr, ) 


-YK^IwKp" exp(/ryr;/,"x,) + 'ZmlC"j„A"p'; exp{i(vm"x,\ 


(2.44) 


<T«I 


(r-4 


where C'ikjm and C"ikini, are the elastic coenicients of the first and second media. The 
boundary conditions must be satisfied not only for all times t, which is a direct result of the 
common term Qxp{io)t) in all equations, but also at all points (xi, X 2 , X 3 = 0 ) on the 

interface. Since the exponential functions Qxp{ia)m1'x,) are linearly independent for a= 0, 

I, . . , 6 Eq, (2.43) and Eq. (2.44) will be satisfied only if all seven exponential functions are 
identically equal. Snell's law follows immediately from this simple fact since this condition 
requires that 


m'lxi - tn',x, 


in]x, = 


m,x, 


mjx, 


nP,x, 


(2.45) 
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I Icic the coordinate vector x lies on the interface plane therefore 

b • X = 0 (2.46) 

where b is normal to the interface as shown in Fig. 4,5. Another way of looking at Eq. (2.45) 
is that for any two different p and y values from the set (0,1, . . ,6) 

(m^ -lll’')-X = 0. (2.47) 

Comparing Eqs. (2.46) and (2.47) shows that the difference m*’ - m’' is always parallel to b, 
i.e., the plane defined by any two of the seven slowness vectors is always perpendicular to 
the interface. By definition, the slowness vector of the incident wave lies in the plane of 
incidence, therefore all the other surface vectors must also lie in the plane of incidence. 
Since the direction of the slowness vector coincides with that of the wave vector, it can be 
conclude that the wave vector of the incident wave as well as all the other wave vectors of 
the reflected and refracted waves lie in the same plane which is normal to the interface', i.e., 
the plane of incidence. Furthermore, 

(in^ - itr)x b = 0 

or , m ^ X b = iir X b, (2.48) 

i.e., all vector products of in'* are equal, which is an alternative form of Snell’s law. This 
relationship can be converted into a more familiar form by introducing polar angles (}»” 
between the interface normal b and the different wave nonnals n™ for the incident wave and 
all refracted and rellected waves so that b. n“ = cos (j)'*. 1 hen, according to the definition of 
the vector product 


in'* X b = /;i'’sin^'* (2.49) 

Since b is a unit vector and a"l, 2..., 6 . 

Fxiuation (2.49) expresses generalized Snell's law for arbitrarily anisotropic media. It is 
important to note that the slowness values m" in general depend on the angles 0™ due to 
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anisolropy, F igure 2.8 illusttnlcs the fuiiclntncntnl principle llutl slowness vector of the 
incident wave and nil the other slowness vectois of the rcllcctcd and refracted waves lie in 
the same plane which is normal to the interface, i.e., the plane of incidence, and that they all 
have the same projection on the plane of the Interface. The same conclusion holds for the 
projections of the wave vectors since k™ = in'*(o. In terms of the wave vectors, we can 
summarize Snell's law as follows: 

y 1 he wave vectors for the incident and rellected/refracted waves all lies in a single 
plane, which also contains the interface normal and is called the incident plane. 

> All their projectiotis on the interface are equal. 



Figwe 2.0 Geneialnrd Snell's law for aibitraiily 
Anij':'hop’i': media 


2.9 CLOSURK 

In this cliaptcr, basic concepts of ray propagation in anisotropic media has been discussed. 
Snell’s law for refraction and reflection phenomena of ultrasonic waves at a planar interface 
has been discussed. 
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CHAPTER 3 
RAY TRACING 

3.1 INTRODUCTION 

Ray tracing is the most important part of a tomographic procedure. The accuracy and 
convergence of image reconstruction is solely dependent on resolution of the ray tracing 
procedure. In the present work, energy is assumed to behave as bundles of rays, each ray 
acting independently from its neighbors. As mentioned in previous Chapters, for 
anisotropic material, direction of energy propagation associated with each ray is generally 
not parallel to the phase velocity direction, the phase velocity defines the orientation of the 
phasefiront associated with each ray. Furthermore, if a material is inhomogeneous then a ray 
path is likely to bend because of the continuous refraction as the ray experiences changing 
material properties. 

For tomographic reconstruction, the information of the ray path followed by a ray between 
two specified positions is necessary. The two points generally represent transmitter and 
receiver position. The method of tracing the ray path between these two positions is known 
as two point ray tracing method. 


Transmitter 



Figure 3. 1 An example showing different iterations in a two point ray tracing method. 

For example, if a defect is situated in a inhomogeneous material as shown in fig. 3.1, then 
determination of the ray path that reaches the defect edge will require two point ray tracing. 
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the defect edge, because of refraction at interfaces and bending in inhomogeneous material. 
As shown in Fig.3. 1, a 60° (angle measured from surface normal) ray aimed at the top of a 
defect is skewed, refracted and bent away from the required endpoint. So for determining 
the route by which, energy will reach the defect requires the solution of a non-linear 
boundary value problem, (the boundary values being the required initial and final ray 
position, the non-linearity arising from the distortion of ray path within the material). As 
there is no general analytical technique exists for solving this problem. Numerical iterative 
techniques have to be employed to find an adequate solution to this problem. Although the 
use of iterative techniques are computationally intensive, they are generally versatile in the 
sense that arbitrary inhomogeneity may be considered and as many arbitrary combination 
of different materials can be taken. 

3.2 THEORY 

The basic concepts and formulation, included in this section are taken from [9]. In the 
absence of body forces, the equation of motion for the displacement vector u in a 
homogeneous anisotropic solid reads as follows; 

(V.C.V).u(R.t) - P - 1^ u(R,t) = 0 (3.1) 

ot 

and is explicitly dependent on space R and time t. Here V is the gradient vector, p is the 
material density. 

The elastic stiffness tensor C has the following structure: " ^ 


f A|| + 2/r,|| i> 

l> A| + 
II A, 


I* 

A. 

A t + 2ft » 


(I 0 ') 

Cl 0 0 

(I 0 0 


U Q 0 \ 

1! » (t 

xi 0 

A , 0 Cf 


{3.2) 
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If the grain direction is parallel to the x-direction, the elastic constants are given as 
C 1 1= A.|| + 2p|| , C22=C.^3= A,1 + 2 Pj. , C44= Pjl, C66=P||, Ci3=V. 



Figures 2 Fiber orientation , phase ve locity and group velocity direction 
relative to cartesian coordinate system 

The plane harmonic wave solutions of Eq. (3.1) are; 

= (3.3) 

where to is the circular frequency, are the polarization vectors, m^=m/V^ is the 

slowness vector that shows the direction of phase propagation m with a value inverse to 
the value of the phase velocity Va. The abbreviations a=SH, qSV, qL denote SH-shear 
horizontal wave, qSV-quasi shear vertical wave, qL-quasi longitudinal wave. Solving the 
corresponding linear algebra problem leads to two dispersion relations; 

= P’ ( 3 - 4 ) 

= p(B ± - 4Ay ^)K2A), (3.5) 
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where, 


^ ) + (a • K)^ (Aj. + 2;Uj.)- v{v 

+ yh)) + i^-ky[{A,+2^J{^^-{A,+2Mr)) 

+ iy + M,.n 

B — + A^ + 2fi^ + (a • K)^ (/l^ + 2//^ — {Aj. + 2fdj . )). 

The direction of the energy transport or the group velocity direction is 

(Y,),=KK + F.P'.(a-K)(l-KK).a 




y“r +(>“/. -/“rXa-K)' 



1 

2 


/ 

+2(a-K)'z^ ± 


F^zSgMfL] 

-JB’-4A , 


yA = +^Mr)- Hy + 2//J 

Ml. )iMr. ~ (^ + v“r )) + + Ml Y 

+ ^Mr)][Mi - i\ + ^Mr)] - 
i^ + MiYiMi+^+^Mr) 


(3.6) 

(3.7) 

(3.8) 

(3.9) 


(3.10) 

(3.11) 


(3.12) 


F^ — F^ + [fJj^ — (Aj. + 2fij.')(Aj. + 2//j.)][//jj (Aj^ + 2//j,)] 

+ MiY (Mz +Ai+ 2/i^) (3. 13) 

As obtained from Eq. (3.8), the direction of the group velocity differs from the phase 
direction wi by an additional vector that corresponds to the projection a onto the plane 
perpendicular to m . 
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The sign ‘ + ’ in Eqs. (3.5) and (3.10) corresponds to the quasi-vertical polarized shear 
wave qSV, and ‘ ’ to the quasi-longitudinal wave qL, respectively. Because of 

symmetric configurations it is sufficient to consider the wave propagation in the j:z-plane. 
According to Snell’s Law, all projections of slowness vector ni onto the interface must be 
equal 

=ml\ (3.14) 

and thei r y-com ponents absent 

(3.15) 

Under these conditions, the mx*^” components can be determined from the dispersion 


relations Eqs. (3.4) and (3.5) or the reflected/refracted SH-wave, we have 


= {-bs ± ^{bl - a,c, ))la. 

(3.16) 


(3.17) 


(3.18) 

~ \./^r ~ P 

(3.19) 

For qS and qL waves, the =S,L) components of the slowness vector are 

defined as the roots of a fourth order polynomial 


{mf y (a, + by^ + ) + {mf f • la^a^m^b, + Ic^a] ) 

+ + b^{a] +a!) + 6c,a>;] + < + ^,^1) 

+ (/«f ) • 2a,ajn^{ml{b^ + Ic^a])-^ e^) 

+ ml {a^ + b^al + c/, ) + m] {d^ + e/, ) + /, = 0 

(3.20) 

where 


~ Pi. ^Pt ) 

(3.21) 

b, = \ + ^Pr ) - V( V + IfX , ) 

(3.22) 

= (A + Pi ypi ~ ^pt )) + (’^ + pl 

(3.23) 

+ +2//r) 

(3.24) 
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= -p{K + - (^ + '^Pt)) 


(3.25) 

(3.26) 


To separate the roots corresponding to qSV-wave and qL-wave, the dispersion relation 
Eq. (5) is transformed with two conditions 


P 


m^s - 5 > 0 


'M 

P 


m 


qL 


B<0 


(3.27) 


(3.28) 


The corresponding group velocity vector for the refracted and reflected waves has to 
point from the interface into the respective material. If a roofs pair from Eq. (3.27) or Eq. 
(3.28) is complex, a surface wave will occur and m 2 ^“ =0. 


3.3 RAY TRACING ALGORITHM 

The ray tracing algorithm follows Ogilivy’s [9], with the exception for the selection of 
the boundary for the ray refraction at ray propagation interval (step 4). The procedure has 
to be followed as described. 

1. For the given Transmitter and Receiver positions, find the group velocity angle, 
for the ray, which will reach Receiver position directly if specimen has no defects. 

2. Find the phase angle using equations 3.8-3.13, for the group velocity angle, using 
iterative procedure, find the group velocity direction and the grain orientation. 

3. Step one time interval forward in the group velocity direction and check that no 
region change has occurred. 

4. If the region changes, reduce the time interval to reach the interface between 
regions and continue with step 5. 

5. In the interface, determine the defect boundary, approximate it as a curve. 
Calculate the new phase velocity directions using equations Eqs. (3. 15)-(3. 16), 
(3.20). 
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6. Return to step 2 and repeat until the specimen boundary is reached. 

7. It distance between final position and receiver position is less than considered 
accuracy-, then terminate, else continue with step 1. 

3.4 EDGE DETECTION 

As the ray reaches interface boundary, the information about the defect boundary is 
required. According to Snell s law, projection of slowness vector on to interface must be 
equal. To identity the interface boundary, all the neighbor pixels to the interface pixel, 
are taken, and then pixels present in interface boundary are identified. Interface boundary 
is identified by material properties of the pixels. A curve, passing through these pixels, is 
approximated.. Then equation of tangent to the curve is calculated. All projections are 
considered along the tangent. Figure 3.3 shows tangent to the curve at the interface pixel 
on interface boundary. 

Figure 3.4 shows magnified view of the interface boundary. Pixels at interface are shown 
by blue line .Red line shows the approximated curve. 
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ray path 

parent material 
defect material 
tanget 

— approximated curve 
pixels at interface 


3.5 CLOSURE 


In this Chapter, basic concepts related to ray tracing are discussed. Then basic steps 
involved in the algorithm have been discussed. The process of edge detection is 
described. 
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CHAPTER 4 

TOMOGRAPHY 


4-1 INTRODUCTION 

Tomography is a method of reconstructing a function from a set of its projections. These 
projections are line integrals of the unknown function computed along certain prescribed 
directions. Data in the form of line integrals arise naturally in measurements involving 
radiation sources as probes. 

Tomographic algorithms are broadly classified as 

1 . I ransform methods — these are direct and hence computationally efficient, but 
require a large number of projections to recover the solution with some degree of 
accuracy. 

2. Series expansion method - these can produce meaningful answers with a fewer 
number of projections, but they are iterative and computationally intensive. 

3. Optimization methods - the optimization approach requires the selection of a 
suitable functionals that has to be extremized, under the constraint of the 
experimentally recorded projections. The performance of this method in terms of 
accuracy and speed depends on quantity to be optimized. 

4.2 CONVENTIONAL ART ALGORITHM 

In tomograpliy, projection data (p) represents an integrated information resulting from 
interaction of source radiation with the material, and is expressed as [1 1] 

P = (4-0 

where and pi are usually taken to be the slowness field in the object and time-of- 
flight (l OF) respectively. Here s represents the ray path of the ray. 

The iterative reconstruction methods of tomography are based on the discretization of the 
cross sectional plane by a square grid. The length of the intercept of an ray with they^* 
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pixel in a given projection is called th weight fraction ( w,) .If is fte field value in the 
f pixel, the f projection is given as 

N 

^Ufj (4 2) 

M ^ ^ 

where M is the total number of rays and N is the number of pixels. 



Figure 4.1 Weight function w-^ for i*^ ray and j* pixel. 

It is required to collect the projection data over a number of view angles to carry out 
tomographic reconstruction. When ultrasound wave is used as a source radiation, f(x,y) 
and px are usually taken to be the slowness field (inverse of sound velocity) in the object 
and time-of-flight (TOF) respectively. The aim in tomography is to obtain f(x,y) from 
given pi. The conventional algebraic reconstruction technique (ART) involves four major 
steps: 

1 . initial assumption of object field 

2. calculation of correction 

3. application of the correction 

4. test for convergence 

The word “convergence” is used here as a stopping criterion. It is not used in the strict 
mathematical sense, where a formal proof is required to show that the solution obtained 
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numerically is lire real solution, rypically, in Gordon ART, the field values jt; in pixels 
are updated for k"' iteration using correction term as 



( Pi -< w\ f*-’>)/(w\ w‘) 


(4.3) 


where < > represents dot product. A, is relaxation parameter and w' is a vector of 
intercepts or weights along ray. 

In isotropic materials, the field value// in a pixel gets updated from all the rays passing 
through that pixel and thus, each has a unique fj , irrespective of direction of ray. 

In case of anisotropic materials, the field value fj is directionally dependent, for example, 
sound velocity in composite materials. So, the value fj in one direction cannot be 
modified by the correction term of ray in another direction. Hence, the f\Q\dif(x,y) is 
expressed as a function of direction ,6 and a set of independent parameter X. it is 
analogous to representing the velocity in composite materials as a function of their elastic 
constants arid the direction of ray. Also Eq(4.2) is modified for anisotropic materials to , 

(4-4) 


where is the field value in /* pixel along i'* ray ,for direction 0i and xj is a vector 
of independent parameter for pixel j. These projection data pi are a result of interaction 
of component of the field function f(x,y) with the source radiation along the direction 0i. 
flius, the problem of tomographic reconstruction for anisotropic materials reduces to the 
tecunsttuction of X in place of f(x,y). 

4.3 MART AIXiORirilMS 

'fhe dilTerence between ART and MART algorithms lies in the manner in which the 
corrections are applied to guessed field values. The corrections are additive in former, 
while multiplicative in the latter. These algorithms are essentially iterative in nature, but 
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inlermediale steps may also involve repetitive calculations in the form of FOR loops. To 
identify the iterative loop, sterL and close.labels with statement numbers have been 
indicated in the description of each algorithm. 


The MAR f algorithms are first summarized. Let (j)ig be the projection due to ray in 0 

direction of projection and / ; be the initial guess of the field value. The projection } 
using the current field value can be obtained numerically [20] as : 

Z l,2,...,Me (4.5) 

y.i 

where Mo represents number of rays along the 6 direction. The parameter to be varied in 
each loop is indicated in brackets. 

Start: 1 vStart iterations (k): ' 

St art: 2 For each projection angle (0): 
vStart:3 For each ray (i0): 

Compute the numerical projection {Eq(4.4)} 

Calculate the correction as : 

A<t>ie = (|)io / ^io 

Start:4 For each cell ( j ) 

If H’lo j is nonzero then : 

MARTI: 

(4.6) 

MART2: 

fT-ffxl (4.7) 


MARO: 

(4.8) 

where p is relaxation factor. 

close:4 

close:3 
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olose:2 

Check (or convergence as: 

If 

abs((/*''-/;/(/^^')] X 100<e 
where c is a suitable stopping criterion. STOP: 
lilsc : Continue 
Close: 1 

Step 3 and 4 (otm the essence of the reconstruction algorithm. Alt three versions include 
the relaxation (actor p. 1 ypical values for relaxation factor reported are in the range 0.01- 
0 . 1 , larger values leading to divergence. It is to be noted that correction calculated in step 
3 is the ratio of the recorded projection data ((j)io) and that calculated from guessed field 

namelyd ^ ) ,which is being iterated. The three versions of MART differ in the manner 

in which the coircctions arc implemented. In MARTI, the weight function is prescribed 
in bitiaty form, being unity if a pixel ray passes through a pixel and zero otherwise. In 
MAR'r2 and MAR 13, the weight function is calculated precisely as the ratio of the 
length of the ray ifiterccpted by the pixel and the maximum dimension of pixel enclosed 
in it 

4.4 rOMOCJRAPH Y AS APPLIED TO COMPOSITE MATERIALS 
Tomographic image reconstruction process can be applied to determine various elastic 
constants of orthotropic composite materials. The process is better known as vector 
tomography, in vector tomography, more than one parameters, in vector form are 
considered for tomographic reconstruction. These parameters are mainly elastic 
constants, density and fiber orientation. Minimum number of projections through each 
pixel should be ct)ual or more than parameters considered for reconstruction. Through the 
di (Terence iti time of fliglit data in actual and estimated ray path, correction is applied on 
group vekKily for all view angles considered. A system of nonlinear equation is obtained 
in terms of group velocity and elastic material properties. Corrected material properties 
are dctermiited by solving the system of nonlinear equations. The process is repeated(^^ 7 
till convergence criteria is satisfied. The method is used for determination of principal 
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elastic constants, mainly C„, and density. The algorithm for the process is as 

follows: 

Step I. Initial assumption of object field 

Input lime of flight data and initial guess for elastic constants. Decide 
convergence criteria. 

Step 2. C'alculation of correction 
Start 1 for each view angle 

■ Start! For each ray ’ 

■ Assume view angle for the ray as group velocity angle. 

■ hind phase angle for given group velocity angle using equation 3.8. 
Ferform two-point ray tracing for given transmitter position and 
receiver position 

• Find intercepts for each pixel, present in the ray path. 

• Find group velocity correction for each pixel, using ART algorithm(eq 
4.3). 

Clo se 2 
CloscJ. 

Step 3 Application of correction 
Start.! For each pixel 

■ Solve for elastic constants, using least square nonlinear method. 

Clos e 1 

Step 4 Test for convergence 

•• If average error is less than convergence criteria or iteration is less than 

1 

maximum number of iterations continue with step 1 else stop. 

In step 3, to solve a system of nonlinear equations, MATLAB inbuilt function ‘fsolve’, is 
used . The program is written in MATLAB having 16 function files and around 2400 
lines of executable statements. A detailed documentation of the program and various 
subroutines is given in Appendix A. In the next section a detailed flow chart is presented 
for the present work. 
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4.5 fLOW DIAGRAM FOR lOMOGRAPHIC IMAGE RECONSTRUCTION 
PROGRAM 
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4.6< I.OStRE 


In this chapter various tomography algorithms are discussed. Then process of elastic 
constant determination is discussed and a detailed flow chart for the developed model is 

included. 
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Chapter 5 
Results and Discussions 

In this chapter, (he results of tomographic recotistiuction in inhomogeneous and 
anisotropic materials are presented. The results consist of leconstraction images 
devcloi>ed by the present model for different materials. 

The present two-ditncnsional model is applied to identify the parent and defect 
materials in a given specimen using time of flight data. The tomographic algorithm is 
also applied to experimental time of flight data, to determine the material properties of 
given specimen. 1 he input parameters for image reconstruction are time of flight Hata 
attd initial gucs.s. Relaxation parameter considered for all reconstruction is 1.0. 

5.1 Ki:< ONS I RK 1 ION USING SIMULATED DATA 

I line ol flight data consists of transmitter position, receiver position and 
conesponding lime of flight. Pixel size considered for tomography is 1x1 mm and for 

ray tracing it is 0. 1 xO, 1 mm. Edge detection is applied to find the interface when ray 
goes from one material to another material. 

I’timograpluc reconstruction for rectangular, elliptical and circular defect is presented 
using simulated data. To obtain simulated time of flight data, correct material 
properties arc assijmed to each pixel and then ray tracing is performed for different 
view angles, Ilie calculated time of fli^t data is used as an input to tomographic 
reconstniction algorithm. Figure 5. 1 shows the ray tracing performed on the specimen 
with rectangultir defect for one of its view angles considered. Material properties for 
parent a^id defect material arc given in section 5.1.1. Transmitter positions for the 
rays arc along z-axis and x-axis and corresponding receiver positions are along z =200 
and X 200 lines, 1116 bending of ray can be observed, as it inters into defect. The 
change in group velocity angle is dependent on properties of parent and defect 
material, incident angle and defect boundary. 

5.1.1 Rectangular defect 

Toiimgraphic reconstruction is applied to construct the images of a specimen with 
rwtangular defect. The size of the ^ecimen is 20x20 mm. The parent material is 
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specimen with rectangular defect 



0 50 100 150 200 

X — > 


Figure 5. 1 Ray paths for group velocity angle 16.7 degree 


Graphite/Epoxy with material properties Cii=145 GPa, Ci2=10.20 GPa, 
C22=C33=13.50 GPa, C44=3.40 GPa, C66=6.80 GPa, density = 1600 kg/m^ and fibers 
are oriented along x-axis. The size of the defect is 6x4 mm and defect properties are 
taken to be Cii~220 GPa, Ci2~10-20 GPa, C22~C33=7.0 GPa, C44=3.40 GPa, C66=6.80 
GPa, density = 2500 kg/m'^ and fibers are oriented along x-axis. Initial guess is taken 
as Cu=120 GPa, Ci2=10.20 GPa, C22=C33=10 GPa, €44=3.40 GPa, C66=6.80 GPa, 
density = 1500 kg/m . Parameters to be reconstructed are Cn, C33 and density. 
Rectangular defect used for simulation is shown in figure 5.1 and 5.2. View angles 
used for simulation are 0°, 90°, 16.7° and 73.3°. Initially it was desired to take view 
angles as 0°, 90°, 15° and 75°, but to receive rays at the center of pixel on the 
specimen boundary; the view angles are slightly modified. The number of projections 
on a pixel is more than parameters considered. 

Figure 5. 3-5.5 shows final reconstructed image after 5 iterations The average error, 
after each iteration is shown in figure 5.6. Average error is calculated as 


Average error = 


/•estimated tactual 

1 ^ A ~Ji 


xlOO 


(5.1) 
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where N==total number of pixels, 

jestmated_ ggjgjjiated material property for i'*' pixel after the 
iteration 

j-actmi _ material property for i*^ pixel 

The maximum error for Cn, C 33 and density are 24.14 percent, 65.8 percent and 14.6 
percent respectively. The maximum error in C 33 is very high, as some of the defect 
pixels are not identified clearly. The defect is approximately identified by above 
images. From reconstructed images it can be observed that a ray passing through the 
defect influences the pixels of parent material, present in its ray path. Also, the 
convergence is observed much faster in parent material than defect material. Effect of 
view angles on reconstructed image can be observed easily. It can be seen that 
average error decreases with successive iterations. Density shows a faster 
convergence as compared to other elastic constants. 


Specimen with rectangular defect 
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Figure 5.3 Reconstruction based on Ci i for rectangular defect 



X — > 

Figure 5.4 Reconstruction based on C 33 for rectangular defect 
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Figure 5.5 Reconstructed based on density for rectangular defect 



Figure 5.6 Variation of average error with iteration 
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5.1.2 Elliptical defect 


Tomographic reconstruction is applied to the specimen with elliptical defect 
(figures.?). The size of the specimen is 20x20 mm. The x-axis and y-axis represents 
pixel number, firom origin, along x-axis and z-axis respectively. The center of the 
ellipse is at (10, 12) with major radius of 5 mm and minor radius of 3 mm. The parent 
material is Austenit with material properties Cii=220 GPa, Ci2=145 GPa, 
C22=C33=150 GPa, C44=i29 GPa, C66=129 GPa, density = 8100 kg/m'^ and grain 
orientation is along x-axis. Defect properties are taken as Cii=80 GPa, Ci2=145 GPa, 
C22=C33=70GPa, € 44=129 GPa, C66=129 GPa, density = 6500 kg/m'’, grain orientation 
along x-axis. Initial guess is taken as Cn=175 GPa, Ci2=145 GPa, € 22 = 033=110 GPa, 
€ 44=129 GPa, €66=129 GPa, density = 7500 kg/m'’ and grain orientation along x-axis. 
Parameters to be reconstructed are €n, €33 and density. It is proposed to have view 
angles of 15°, 75°, 105° and 165° but to make the rays reach at the center of pixels the 
actual angles are taken to be 16.7°, 73.3°, 163.3° and 106.7°. Figure 5.8 shows ray 
path for view angle 73.3°, used for simulation. Edge detection is applied after three 
iterations. Effect of edge detection can be observed m Figure 5 . 8 . 

The final reconstructed images for € 11 , €33 and density are shown in figure 5.9, 5.10 
and 5.11 respectively. Figure 5.13 shows the average error for all these properties .It 
can be observed that, in this case, density converges faster as compared to other 
properties. Figure 5.12 shows the reconstruction of €n at various stages. The shape 
of the figure can be identified in 2 ““^ iteration itself but the error in €n is high and a 
drop in error iBrom 1^ iteration to 6 **^ iteration can be observed. After 6 *^ iteration the 
€11 values for parent as well as defect pixels converge to actual values. Effect of view 
angles can be figured out from these reconstructed images. It can be observed that, 
after 3 "^ iteration, application of edge detection, results into faster convergence. 

The maximum errors, after 6 *^ iteration, for €u, C 33 ands density are 27.3 percent, 
39.8 percent and 24.61percent respectively. Error in €33 is high, as some defect pixels 
are not identified exactly. 
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Figure 5.9 Reconstruction based on Cu for elliptical defect. 


Reconstruction based on C33(Pa) ^ qIo 
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X — > 

Figure 5.10 Reconstruction based on C 33 for elliptical defect. 
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Figure 5.11 Reconstruction based on density for elliptical defect 


Reconstruction based on C1 1 (Pa) after 1 st iteration q 1 i 
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Figure 5.12(a) Reconstruction based on Cn after 1®* iteration for elliptical defect 
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Figure 5. 12(c) Reconstruction based on Cn after 3”* iteration for elliptical defect 
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Reconstruction based on C1l(Pa) after 4th iteration 
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Figure 5. 12(d) Reconstruction based on Cn after 4 iteration for elliptical defect 
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Figure 5.12(e) Reconstruction based on Cn after 5*^ iteration for elliptical defect 
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Figure 5.13 Variation of average error with iteration 
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5.2 RECONSTRUCTION USING EXPERIMENTAL DATA 


The performance of the proposed tomographic algorithms in the context of 
experimental data is reported in the present section. An experimental set up has been 
developed in the laboratory by Mr. Gola Dalliraju and Mr Sarin, whose details are 
given in appendix B. The specimen size is 30x30 mm, with square defect of size 
lOx 10 mm in center (Figure 5. 14). The parent material is glass/epoxy and the defect is 
of Aluminum. For the considered specimen, volume fraction for the glass fiber is 0.45 
(determined by Bum testing performed by Yekambar Rao, in Fracture Mechanics 
Lab). For the volume fraction 0.45, the material properties for Glass/epoxy are found 
to be Cn=39.17 GPa, € 33 = 8.4 GPa, Ci3=2.18 GPa, C66=4.14 GPa and density =1800 
kg/m^. For defect, grain orientation is considered along y-axis. The experimental data 
is used for image reconstmction based on Cu only. The view angles, 0° and 90° are 
used for experiments. Initial guess considered for reconstmction is Cn=35 GPa. 
Experimental data used for reconstmction is included in Appendix B. 

For tomographic reconstmction time of flight data is determined for both view angles. 
The reconstmcted image after 8 iterations is shown in figure 5.15. Average error for 
defect, considering actual values as, Cii=72.5 GPa is determined. Average error in 
determination of Cn is found to be 12.14 percent. After 8 iterations, the average 
value of Cii for parent and defect material is found to be 44.6 GPa and 63.7 GPa 
respectively. 

Reconstmcted images for experimental data showed similar kind of convergence as 
shown by simulation. However the error involved in image reconstmcted from 
experimental data is higher. Errors in experimental data can be associated with 
misalignment of the laser beam and the inaccuracy in image processing operations. 
The reconstmction quality obtained from experimental data is satisfactory. The top- 
right edge of defect is not detected properly by the algorithm, due to presence of noise 
in experimental data. Reconstmction image for Cu (Figure 5.15) shows that effect of 
rays at 0° is dominating. Pixels present in right side of defect do not show proper 
convergence for reconstmction of parameter Cu due to noise in experimental data. 
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Figure 5.14 


Figure 5.15 Image reconst 


5.3 EFFECT OF VARIOUS PARAMETERS ON THE PERFORMANCE OF 
ALGORITHM 


In this section, performance of the algorithm as a function of different parameters is 
studied. 

• Effect of initial guess 

Initial guess of the elastic constants have a significant effect on theconvergence of the 
parameters (table 5.1). Convergence is poor if difference in initial guess and actual 
value is high as shown below. 


Table 5. 1 Effect of initial guess on average error (after 5 iterations) 


Initial guess/actual value (%) 

10 

50 

80 

90 

110 

150 

Average Error (%) 

28.2 

16.7 

4.7 

2.5 

3.1 

15.9 


• Effect of relaxation parameter 

The relaxation parameter influences the rate of convergence. Figure 5.16 shows 
average error observed for different relaxation parameter. There is no significant 
difference in the quality of results for all values except at very low value. The value of 
relaxation parameter recommended is between 0.8 and 1.2. 

• Effect of number of Projections 

As the number of projections through a pixel is taken more than parameters for 
reconstruction, a faster convergence is observed. Figure 5.17 shows how average error 
reduces when the number of projections taken in a pixel is increased. It can be is 
observed that rate of convergence is faster for more number of projections but after 
certain limit, no significant improvement could be observed. Generally, for n 
parameters, 2 n projections through a pixel give optimum results. 

• Effect of edge detection 

On application of edge detection results show better convergence. Figure 5.18 and 
5.19 shows reconstmcted images for Cn without applying and applying edge 
detection. The average error without considering edge detection is 8.6 percent, after 
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5* iteration. While on applying edge detection, the average error reduced to 3.8 
percent. 




Figure 5.17 Effect of Number of projections throu^ a pixel on performance of the 

algorithm 
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Reconstruction based on C33(Pa) without considering edge detection 


Figure 5.18 Reconstructed image for C 33 without applying edge detection 


Reconstruction based on C33(Pa) after applying edge detection 
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Figure 5.19 Reconstructed image for C33 after applying edge detection. 



5.4 CLOSURE 

In this chapter, present model is applied on simulated and experimental data and 
results are summarized. Effect of various parameters on the performance of 
reconstruction algorithm is discussed. 
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Chapter 6 

Conclusions and Scope for Future work 

6.1 CONCLUSIONS 

A tomographic image reconstruction model, applying ray tracing, has been developed. It is 
applied to determine principal elastic constants and density of composite materials. The 
model is very useful to predict and understand the effects of various parameters on 
tomographic image reconstruction. From this study, following conclusions can be drawn. 

r- The performance of tomographic image reconstruction model is improved when 
number of projections considered from a pixel is much more than number of 
parameters considered for optimization. 

y Application of edge detection on the ray tracing algorithm, improves the resolution of 
reconstructed image and average error is reduced. 

> It has been found that if a specimen is scanned across a surface at constant beam 
angle, then there is a possibility that ultrasonic ray energy may not penetrate certain 
regions. Therefore in such cases to ensure that the whole component can be 
examined, more than one probe angle may be needed. 

^ Application of ray tracing in tomographic reconstruction, insures a goo4 rate of 
convergence. The model is found to be successful in the reconstruction of images and 
determination of material properties, for rectangular as well as curved defects. 

6.2 SCOPE FOR FUTURE WORK 

> The method can be extended for determination of shear modulus, rigidity modulus 
and fiber orientation of composite materials. 

> The model can be extended for a combination of quasi-longitudinal and quasi- 
shear waves for identification of defects. 

> With few modifications, triangular pixels can be taken, instead rectangular. This 
will improve the accuracy when defects are in curved shape. 
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Appendix-A 


This Appendix will describe the file structure of the program developed in the present 
work and it will discuss different fimctions used in the program. 

File structure 

In whole 16 function files are used in the program. Each fimction file performs one or 
more specific tasks. The basic structure of the whole program based on these files is 
given in the form of flowchart below. 



Functions used in the program 

mainjomography.m : It is the main program. Different functions are called in it and it 
determines the elastic constants for each pixel after each iteration and updates the pixel 
properties. 
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inputdata.m : It specifies the specimen geometry, no. of pixels, initial guess for elastic 
constants. Time of flight data including Scanner position. Receiver position are also 
specified in it. In addition, it specifies step size for ray tracing. 

find_weightage.m : for a given Scanner position and Receiver position, it determines the 
ray which is reaching to receiver position through an iterative procedure. Then it calls 
function find_intercept.m to calculate the intercept in each pixel for the ray. 

find_finaljpixel.m ; It determines the final pixel of a ray, for given Scanner position and 
group velocity angle. It stores the ray path and pixel properties for each pixel in the ray. 

check_boundary.m : It checks if ray has reached the specimen boundary or not. 

function_phase_angle: For given material properties and group velocity angle it 
determines the phase angle. 

function_tangent_direction.m : For given interface pixel, it approximates a curve and find 
the tangent and normal direction at the interface. 

slo'wness_vector.m; It calculates the slowness vector for given material properties, group 
velocity angle and phase angle. 

slowness_vector_exit.m: It calculates the slowness vector, when ray reaches an 
interface,and refracts according to change in material properties. 

find_intercept.m : For a given ray path, it determines the intercept length for each pixel, 
and calculates the corrected group velocity using ART algorithm. It stores the 
information required for elastic constant determination for each pixel. 
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group_vel_direction.m: It calculates the group velocity direction for given phase angle 
and material properties. 

roots.m ; It calculates the roots of 4 degree equation required for finding the slowness 
vector at interface. 

find_elastic_constants.m : For a given cell, it determines the corrected elastic constants 
from the corrected group velocity data. It uses Least square nonlinear Matlab inbuilt 
function to solve for elastic constants. 

GV.m : It calculates the error in corrected group velocity and present group velocity and 
supplies it to find_elastic_constant.m. 

error_estimation.m: It calculates the percentage average error in elastic constants at 
respective pixel , after each iteration . 

tomography.m It calculates the time of flight data for given Scanner position and 
Receiver position , which is further used as simulated data for image reconstruction. 
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Appendix-B 
Experimental Set Up 


This Chapter briefly describes the details of laser based ultrasonic setup and data 
collection. 

EXPERIMENTAL SETUP: LASER BASED ULTRASONICS 

Nd: YAG Pulse Laser is used to generate ultrasonic waves in any media, which is to 
be inspected. Heterodyne type Laser Interferometer is used to detect the transmitted 
wave, the signals are then amplified and digitized using a Yokogawa DL1740 digital 
oscilloscope. The oscilloscope is triggered using a synchronization signal from the 
pockels cell of the pulsed laser. Recorded waveforms are transferred to a HP Pentium 
III computer over a USB/Ethemet interface for subsequent storage and analysis which 
is then captured in a Digital Storage Oscilloscope. The schematic layout of 
experimental setup is shown in Fig.B.l and a photograph of the same is shown as Fig. 
B.2. Small part of the beam is allowed to fall on the specimen to study the Nd: YAG 
laser beam characterizations. For this purpose the Teflon material with different slits 
viz. rectangular of 10mm x 1mm, circular slits of 1 1mm, 1 2mm, 1 3 mm and 
circumferential slits of 1mm width at the mean radii at 4mm, 5mm, 6mm were used. 
All these slits are made in 1 25mm Teflon material. A brief description of the setup is 
given in the following paragraphs. 

<D Nd: YAG Pulsed Laser Ultrasonic Generator 
The 5000 DNS series pulsed Nd: YAG laser is built on a modular concept. The 
optical configuration allows variable setting of the optical parameters. It consists of 3 
major components, the laser-head, the power supply and the cooling unit. The heart of 
the system is the pumping structure, which houses the Nd: YAG rod and the flash 
lamp. The lasers are built on an electro-optically Q-switched oscillator. This oscillator 
uses a pockels cell Q-switch to produce pulses of high intensity and short duration (5- 
7 ns). Q- Switched Lasers are often used as a non-contact ultrasound source in non- 
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destructive testing of materials. Q-switched lasers typically have nanosecond (ns) 
pulse durations and generate broadband ultrasound waves, though longer laser pulses, 
of 100 microseconds or greater, have been used for NDE. A variable reflectivity 
output coupler allows the extraction of high energy on a single spatial transverse 
mode. This leads to laser beams of low divergence, and to high conversion 
efiSciencies in the harmonic wavelengths (1064, 532, 355, 266 nm). Different 
harmonic generators extend the wavelength range to the second, third and fourth 
harmonics. 

The active medimn of the laser operates on transitions of triply ionized 
Neodymium atoms (Nd^"*^), which take place of another ion (Yttrium) in the host; 
Yttrimn Aluminum Garnet known by the acronym YAG. The laser operates as a 4- 
level system. The intense broad-spectrum light of the flash lamp populates the upper 
level. Once in the higher energy level, the Neodymium ions drop to a metastable 
level, producing a population inversion. The lower level decays by a fast non- 
radiative process to the grotmd state. The strongest Neodymium line is 1064 nm. 

■ Optical Heterodyne Laser (He-Ne) Probe 

The SH-130 probe is designed to measure transient mechanical displacements of very 
low amplitude. It is specially devoted to measuring displacements generated by the 
propagation of an acoustic or ultrasonic wave. The system consists of a compact 
optical head and an electronic signal-processing unit. The optical head integrates high 
stability, low power laser somrce fijr fest detection with a high spatial resolution. 

The electronic signal processor delivers a response proportional to the 
displacement of the target, with a high bandwidth. The output signal is automatically 
calibrated to give the absolute value of the measured displacement. The principle of 
detection (heterodyne interferometry) makes the system insensitive to external 
vibrations. The compactness of the system allows for a wide range of operating 
conditions. 

■ Digital Storage Oscilloscope (DSO) 

The setup utilizes a Yokogawa DL1740 (four channel, one GSa/sec, 500 MHz) 
Digital Storage Oscilloscope, having built in Zip drive, Ethernet, USB, GPIB and 
Serial Ports for communication with external PC’s/ systems. 
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■ Computer 

The DSO is interfaced to a Pentium III based PC and Data is stored on-line using 
communication through USB and Ethernet ports. 
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Figure B1 Experimental Set Up 
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Fig.B.2 Photograph of LBU set-up 
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EXPERIMENTAL DATA USED FOR RECONSTRUCTION 


X co-ordinate 
(transmitter) 

y co-ordinate 
(transmitter) 

X co-ordinate 
(receiver) 

y co-ordinate 
(receiver) 

time of flight 

(mm) 

(mm) 

(mm) 

(mm) 

(sec) 

0 

5 

30 

5 

6.78E-06 

0 

15 

30 

15 

6.78E-06 

0 

25 

30 

25 

6.78E-06 

0 

35 

30 

35 

6.78E-06 

0 

45 

30 

45 

6.78E-06 

0 

55 

30 

55 

6.78E-06 

0 

65 

30 

65 

6.78E-06 

0 

75 

30 

75 

6.78E-06 

0 

85 

30 

85 

6.78E-06 

0 

95 

30 

95 

6.78E-06 

0 

105 

30 

105 

6.84E-06 

0 

115 

30 

115 

7.30E-06 

0 

125 

30 

125 

8.1 IE-06 

0 

135 

30 

135 

7.87E-06 

0 

145 

30 

145 

7.70E-06 

0 

155 

30 

155 

7.53E-06 

0 

165 

30 

165 

7.37E-06 

0 

175 

30 

175 

7.21 E-06 

0 

185 

30 

185 

7.04E-06 

0 

195 

30 

195 

7. 11 E-06 

0 

205 

30 

205 

6.78E-06 

0 

215 

30 

215 

6.78E-06 

0 

225 

30 

225 

6.78E-06 

0 

235 

30 

235 

6.78E-06 

0 

245 

30 

245 

6J8E-06 

0 

255 

30 

255 

6.78E-06 

0 

265 

30 

265 

6.78E-06 

0 

275 

30 

275 

6.78E-06 

0 

285 

30 

285 

6.78E-06 

0 

295 

30 

295 

6.78E-06 

5 

0 

5 

30 

1.46E-05 

15 

0 

15 

30 

1.46E-05 

25 

0 

25 

30 

1.46E-05 

35 

0 

35 

30 

1.46E-05 

45 

0 

45 

30 

1.46E-05 

55 

0 

55 

30 

1.46E-05 

65 

0 

65 

30 

1.46E-05 

75 

0 

75 

30 

1.46E-05 

85 

0 

85 

30 

1.46E-05 

95 

0 

95 

30 

1.46E-05 

105 

0 

105 

30 

1.16E-05 

115 

0 

115 

30 

1.16E-05 

125 

0 

125 

30 

1.16E-05 

135 

0 

135 

30 

1.16E-05 
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X co-ordinate 
(transmitter) 

y co-ordinate 
(transmitter) 

X co-ordinate 
(receiver) 

y co-ordinate 
(receiver) 

time of flight 

(mm) 

(mm) 

(mm) 

(mm) 

(sec) 

145 

0 

145 

30 

1.16E-05 

155 

0 

155 

30 

1.16E-05 

165 

0 

165 

30 

1.16E-05 

175 

0 

175 

30 

1.16E-05 

185 

0 

185 

30 

1.16E-05 

195 

0 

195 

30 

1.16E-05 

205 

0 

205 

30 

1.46E-05 

215 

0 

215 

30 

1.46E-05 

225 

0 

225 

30 

1.46E-05 

235 

0 

235 

30 

1.46E-05 

245 

0 

245 

30 

1.46E-05 

255 

0 

255 

30 

1.46E-05 

265 

0 

265 

30 

1.46E-05 

275 

0 

275 

30 

1.46E-05 

285 

0 

285 

30 

1.46E-05 

295 

0 

295 

30 

1.46E-05 
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